Non-Conventional Anderson Localization in Bilayered Structures with Metamaterials 
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We have developed an approach allowing us to resolve the problem of non-conventional Anderson 
localization emerging in bilayered periodic-on-average structures with ahernating layers of right- 
handed and left-handed materials. Recently, it was numerically discovered that in such structures 
with weak fluctuations of refraction indices, the localization length Lioc can be enormously large 
for small wave frequencies uj. Within the fourth order of perturbation theory in disorder, <^ 1, 
we derive the expression for Lioc valid for any oj. In the limit a; — !> one gets a quite specific 
dependence, L;~^ oc a^w^. Our approach allows one to establish the conditions under which this 
effect can be observed. 

PACS numbers: 72.15.Rn, 42.25.Dd, 78.67.Pt 
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Introduction. - The remarkable success in nano- and 
material science has led to a burst of interest to the 
study of wave propagation and electron transport in pe- 
riodic one-dimensional (ID) systems (see, e.g., Refs. [1] 
and references therein). Among others, the systems of a 
particular interest are bilayered .structures in optics and 
electromagnetics, or semiconductor superlattices, as well 
as alternating quantum wells and barriers in electronics 
0. The interest to this subject is due to a possibility to 
create devices with prescribed properties of transmission 
and/or reflection. 

As is shown in recent studies the new perspec- 

tives are related to specific optic properties of meta- 
materials that can be embedded in periodic structures. 
A particular system that has attracted much attention, 
consists of two arrays of a, 6— layers with equal widths, 
da — db, one of which is of right-handed (RH) material 
with positive refraction index, > 0, and the other is of 
left-handed (LH) material with negative index, < 0. 
In such periodic structure with equal = |nf,| the phase 
shift of wave function gained in a— layer is canceled by 
the subsequent shift in next &— layer. Therefore, the to- 
tal phase shift vanishes after passing N unit (a, b) cells, 
that, together with perfect transmission, results in a kind 
of "non-visible" structure. 

However, even a small amount of disorder can destroy 
the above ideal picture due to an emergence of localiza- 
tion [3]. To great surprise, the numerical data 51 ob- 
tained for the model with weakly disordered refraction 
indices and n^, have demonstrated very fast diver- 
gence of localization length, Lioc ~ for w — >■ 0. In 
order to explain such a non-conventional Anderson local- 
ization, in Refs-d-lzj it was suggested that the RH-LH 
structures can not be fully described via the Lyapunov 
exponent A as an inverse localization length. Recently, 
this study was extended to the case of off-axis incidence, 
for which the transport characteristics have been inves- 



tigated in great detail • 

In Ref. another type of disorder was analyzed ac- 
cording to which both widths da and db are slightly fluc- 
tuating. Specifically, the general case of correlated dis- 
order was analytically studied by assuming any kind of 
statistical correlations in widths of a and 6— layers, as 
well as the inter-correlations between the two widths. It 
was shown that for any ratio between da and db the local- 
ization length is governed by the unique expression, no 
matter whether the a, 6— layers are both of normal ma- 
terials (RH-RH array) or with alternating right-handed 
and left-handed materials (RH-LH array). The same 
statement stems from the analysis of RH-RH and RH- 
LH arrays with fluctuating refraction indices [lo| . how- 
ever, only when da db- These results indicate that the 
abnormal localization emerges in a very specific situa- 
tion, and not only due to inclusion of LH-arrays into the 
structure. 



In order to shed light on this problem, in Ref. [ll| the 
RH-LH model with equal widths da = db has been an- 
alyzed with the use of a powerful method developed in 
Ref. [l^ . It was found that the phase of wave, propa- 
gating through the RH-LH array with fluctuating refrac- 
tion indices, is described by a highly non-uniform distri- 
bution. This effect is somewhat similar to that arising in 
standard ti ght- binding models close to band edges (see, 
e.g. review [13'] and references therein). Further analysis 
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of this fact has led to a remarkable result 
Lyapunov exponent vanishes in the quadratic approxi- 
mation in disorder parameter. Nevertheless, the prob- 
lem has remained open since the evaluation of next order 
terms for the Lyapunov exponent was confronted with 
severe technical difficulties. 

In this Letter we suggest the method allowing one to 
resolve the above problem of non-conventional Anderson 
localization. We derive the expression for the Lyapunov 
exponent within the fourth order perturbation theory. 
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which is valid for any value of frequency w and leads to 
the asymptotic expression Ljoc oc a~^uj~^ for small u). 
Our approach elucidates the mechanism that is responsi- 
ble for such an abnormal localization, and allows one to 
specify the conditions under which it emerges. 

Model. - The model describes the propagation of elec- 
tromagnetic waves of frequency lo through an infinite 
array of two alternating a and b layers (slabs), respec- 
tively specified by the dielectric permittivity Safi, mag- 
netic permeability /ia.b, refractive index na,b = i/£a^h/Vbj 
impedance Za,b — fJ'a,b/na,b and wave number ka,b = 
i^na,b/c. We consider two systems: the homogeneous 
stack when both a and b layers are made of right-handed 
optic materials, and mixed stack when a layers contain 
RH-material while b layers are of left-handed material. 
Following Ref. we admit that the compositional dis- 
order is incorporated via dielectric constants Sa.b only, 

na{n) = 1 + ?7a(n), Ma = 1, (la) 
nb(n) =± [1 + 776 (n)], fib = ±l. (lb) 

Here n enumerates the nth unit (a, b) cell. The up- 
per/lower sign stands for the RH/LH material respec- 
tively. Without disorder, T]a^b{n) = 0, the RH-RH struc- 
ture is just the air without any stratification, while the 
RH-LH array represents the ideal mixed stack (eq — 
Mq = 1, efc = /if, = — 1) with perfectly matched slabs 
{Za = Zb = 1). All slabs have equal thicknesses, 
da = db ^ d/2, thus the size of any (a, b) cell is d. 

We specify that delta-correlated sequences r]a,b{n) have 
zero average, {ria,b(n)) = 0, and variance cr^. 



{r]a{n)r]b{n')) = <7^SabS„ 



(2) 



Here (...) stands for the average along arrays, which is 
assumed to be equivalent to ensemble average. Numer- 
ically, to generate rja{n),rjb{n') we use fiat distribution, 
however our expressions are given in general form. 

Basic relations. - As is shown in Refs. 9-11], the 
model can be written in the form of area-preserving map, 

Qn+l — AnQn + B^Pn, Pn+1 — ~CnQn + -D„P„. (3) 

Here the "coordinate" Qn and "momentum" P„ refer to 
the wave function and its normalized derivative, respec- 
tively, taken at left-hand edge of the nth unit (a, b) cell. 
The factors An, i?„, C„, D„ read 

A„ = cos ifa cosipb- Z^^ Zb sin ipa sin ipb , (4a) 
B„ = Za sin ipa cosipb + Zb COS ifa sin ipb, (4b) 
C„ = ^ sin ipa cos ifb + Z^ ^ cos (ySc sin (^Sb , (4c) 

Dn = COS If a COS ifib - ZaZ^^ siu ifia siu (ySfc. (4d) 

For non-zero disorder the coefficients ([4]) depend on the 
cell number (discrete "time") n, due to randomized re- 
fractive indices ([l} entering both the impedances Zafi{n) 



and phase shifts ^Pa.bin), 

ipa{n) ^ ka{n)d/2 ^ ip[l + riain)], (5a) 
(pb{n) = kb{n)d/2 — ±ip[l + rib{n)], (p = ujd/2c. (5b) 

It is useful to pass to polar coordinates Rn and 9n by 
the standard transformation, 

Qn = Rn cos 9n, Pn = Rn sin On. (6) 

The map in the radius-angle presentation gets the form 

(7a) 



(Rn+l/Rn) — dOn+l/dOm 

'Cn + DntaROn 



arctan 



An + Bn tan 6'„ 



(7b) 



The localization length Lioc can be derived via the Lya- 
punov exponent A 11, 1^, 



d 



A 



1/ ddn+l 
'2\^ dOn 



(8) 



Thus, the 0-map (j7bp is the unique equation to be 
treated. We assume that the disorder is weak. 



cr^ < 1 and (aip)"^ < 1, 



(9) 



that allows us to develop a proper perturbation theory. 

Phase distribution. - The phase distribution p{9) can 
be found in the way described, e.g., in Refs. [ij, [l^. By 
expanding the exact 0-map (j7bp up to the second order 
in perturbation and taking into account the uncorrelated 
nature of the disorder, see Eq. ([2]), one can obtain. 



where 



'l-ila{n)U(en)Tr]b{n)U{en 
-a'^Widn), 



7/2) 



(10) 
(11a) 



U(&) — Lp + sin Lp cos(26' — 1^9) , 
W{e) = Lp[cos{2e - 2lp) ± cos(26l - 27)] 
+ sin(p[sin6' sin(0 — ip) ± s\n{9 — 7/2) sin(6' — ip — 7/2)] 
-I- sin^ (^sin(46' — 2(/3 — 7) C0S7 . (Hb) 

The unperturbed Bloch phase shift 7 over a unit (a, &) 
cell is defined as 



7 = (^±(^= (l±l)wd/2c. 



(12) 



where "plus" stands for the RH-RH structure, and "mi- 
nus" refers to the mixed RH-LH array. To proceed fur- 
ther, we pass from Eq. (|10p to stationary Fokker-Plank 
equation for the phase distribution p(&) (l2. 13|. 



{V\Q) + U\d-^I2)\ p(B) 



+2 



d9 



7 



+ w{e) p(0) = o. 



(13) 



which should be complemented by the normalization con- 
dition and the condition of periodicity p{6 -I- tt) = p{9). 



3 




0.8 


c) 


0.8 




d) 


0.6 




0.6 
















q1 0.4 




0.4 






0.2 




0.2 














1 i i 





n/4 nil 3lc/4 t 7t/4 nil 3ll/4 n 

9 e 



FIG. 1. (color online) a) phase space trajectory generated by 
Eq. @ for RH-RH array with N ^ W^, (p ^ 2n/30, for zero 
disorder (solid circle), and for — 0.003 (scattered points); 
b) one trajectory for RH-LH array with A'^ = 10®, (p = 2-k/5, 
cr^ = 0.003. c) p{e) from Eq. Q for RH-RH array (his- 
togram), and Eq. (fTl)) (horizontal hne); d) p{6) from Eq. ((3]) 
for RH-RH array (histogram), and Eq. {TTJ (solid curve). 



One can see that the form of solution of Eq. 
strongly depends on whether the phase shift 7 is non- 
zero (RH-RH array) or vanishes (RH-LH array). 

Homogeneous RH-RH array. - In such a structure the 
Bloch phase (IT^ is finite, 7 = ud/c, and for weak 

disorder the term in Eq. (1131) containing 7/17^ prevails 
over the others for any value of phase shift ip. Therefore, 
the phase distribution is uniform within the first order 
of perturbation theory. 



Pie) - i/tt. 



(14) 



The Lyapunov exponent can be derived by substitution of 
Eq. PH)) into definition ^ , and expanding the logarithm 
up to quadratic terms in disorder. After, the subsequent 
averaging over 9 is trivial and one gets, 

d / Lioc ^ ^ = o''^ sin^ (fi ; ip — ujd/2c. (15) 

This gives standard w— dependence, A oc for w — > 0, 

A = (j^uj^d'^/Ac^. (16) 

Mixed RH-LH array. - The principally different situ- 
ation emerges for the RH-LH array. In this case 7 = 
independently of the phase shift ip. As a result, we have 
W{9) = -U{e)U'{e) in Eq. ([11]), and Eq. ^ leads to a 
highly nonuniform phase distribution. 
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p(0) = -V^2_3inV/{/(0). 



(17) 



Fig. [T] displays perfect agreement between Eqs. (ITi|) . p7)) 
and data obtained by the iteration of the exact map ([3]) . 



The above result means that to calculate the Lyapunov 
exponent via Eq. (|8]), one needs to perform an average 
with the distribution p(9) given by Eq. ([T7|) . Surprisingly, 
the use of Eqs . (jlOp and (fT7|) results in the vanishing 
value of A [ll|. Therefore, the Lyapunov exponent is 
determined by next orders of the perturbation theory. 
However, the direct evaluation of high order terms in 
p{6) is not possible with the above method [l^ . 

In order to proceed further, we suggest another 
method. The numerical data in Fig. [T]d manifest that 
the trajectory has the form of fluctuating ellipse speci- 
fied by angle with respect to axes, and by fixed aspect 
ratio. Therefore, one can expect that in new variables 
Qm Pn obtained by rotating and rescaling the axes Q, P, 
the trajectory transforms into fluctuating circle. Thus, 
the distribution of a new phase in the considered ap- 
proximation will be uniform. 

To follow this recipe, we rotate the Q, P- axes by angle 
T, with further rescaling the axes due to free parameter 
a. In new coordinates the expressions ([3]) and (HI) - (|8]) 
conserve their forms, however, with the factors, 

An = An cos^ r -|- {Bn — Cn) sin T COS T -|- Dn siu^ r, 

(18a) 

Bna~^ = Bn COS^ T — {An — Dn) s'm T COS T Cn Siu^ T, 

(18b) 

Cn<3 — Cn COS^ T + [An — Dn) sin T COS T + Bn SVC? T, 

(18c) 

Dn = Dn COS^ T — {Bn — Cn) Sin T COS T + An Sll? T, 

(18d) 

instead of initial An , Bn , C„ , Dn . 

Now the distribution p{Q) for new phase Q can be 
found starting from the quadratic expansion of Eq. (|7b|) 
with new coefficients ([TSl) and 7 = 0, 



e„+i - e„ = [7?a(n) - m{n)]v{en) + (TV(e„)T/'(e„). 

(19) 

Here the function ^(8) is defined by 

V{Q) = sin (f sin(2r — ip) sin 28 

-f — [<y9 — sin(/5COs(2r — c/?)] [cos 28 — 1] 



1 

'2^ 



[(p -h sin(^cos(2r - (p)][cos28 + 1]. (20) 



The stationary Fokker-Plank equation corresponding to 
8-map (IT9)) reads. 



v\e)—p{e) + v{e)v'{e)p{e) 



0. (21) 



From this equation one gets that the phase distribution 
is uniform, p(8) = I/tt, and the trajectory is, indeed, a 
fluctuating circle provided that 
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V{Q)V'{Q)=Q. 



(22) 
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FIG. 2. (color online) (a) Phase space trajectory in new vari- 
ables {Q,P); (b) distribution p(0) generated by the trans- 
formed map with (|18p and (|23|l . for 7 = 0, ip — 2n/5, 
cr^ = 0.02 and N = lO'^. 



With the use of Eqs. (I^U)) and we now can obtain 
the desired expressions for the angle r, parameter a and 
function V{Q) (which is actually no more 0-dependent), 



2' 



if + sm if 
tp — sin ip ' 



^(©) 



— 1/(^2 



sin^ Lf. 



(23) 

The data presented in Fig. ^ confirm the success of our 
approach: in new variables the trajectory is a fluctuating 
circle and the phase distribution is uniform. 

The Lyapunov exponent A can be now obtained via 
Eq. ([8|) with the change 6n — ©n- Taking into account 
that A vanishes within quadratic approximation in dis- 
order, we expand the 0-map of the form (|7b|) with the 
coefficients (|18p up to the fourth order in perturbation. 
By substituting the resulting expression into Eq. dH) and 
expanding the logarithm within the same approximation, 
after the averaging over 0„ with flat distribution, we ar- 
rive at final expression, 



d 

Lloc 



C + 2 4 [(2(y9^ — sin (/s) cosi^ — (/9sin(/3]^ 
— ; — cr 



9 ■ 2 

— sm ip 



(24) 



Here -2 < ( < 00 

^4\ //„/„\2\2 



stands for the excess kurtosis, 
= (?y(n)'*)/(r7(n)^)^ — 3, the constant specified by the 
form of distribution of rja.b{n). For Gaussian and flat 
distributions we have C = 0, — 6/5, respectively. 
Eq. determines the asymptotics for small <y9. 



d 

Lloc 



= A = 



24 
335^ 



■(C + 2)aV^ for(^<l<f7-\ (25) 



that results in a quite surprising frequency dependence of 
the Lyapunov exponent, A oc cj^. Thus, the dependence 
A (X cj^, numerically found for small uj in Refs. 0, Q 
should be treated as the intermediate one, apparently 
emerging due to not sufficiently large lengths N over 
which the average of A is performed. 

Figure [3] shows excellent agreement between numer- 
ical data obtained for the localization length from ex- 
act Eq. ^ and Eqs. dM]), ^ for RH-LH, as weU as 
Eqs. dlSl), (HH) for RH-RH arrays with C = -6/5. For 
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FIG. 3. (color online) Localization length for RH-LH and 
RH-RH arrays of different lengths N versus the normalized 
wavelength 2-!y/(P — Anc/ujd, for = 0.02. 



sequence lengths N = 10^; 10^ and 10^ the data are ob- 
tained with ensemble averaging over 100 realizations of 
disorder, while at iV = 10^^ only one realization is used. 

Conclusions. - Our approach allows one to resolve 
the problem of non-conventional Anderson localization 
for bilayered RH-LH arrays with equal widths, da = dh, 
and derive Eq. (p4)) for the localization length. As we 
show, the peculiarity of this model is entirely due to zero 
value of unperturbed Bloch phase, 7 = 0. We derive the 
expression for the Lyapunov exponent which appears to 
be defined by the fourth order of perturbation theory in 
disorder, and valid for any value of frequency uj. Our 
results prove that for small w the localization length is 
enormously large, Lioc oc a^'^u^^. 

According to the analysis given above, one can con- 
clude that the non-conventional dependence ijoc oc 
a^'^uj^^ also emerges when da ^ ^6, in case of equal 
unperturbed optic lengths, Uadb = |fib<if)|, provided that 
the impedances are also equal, Za = Zh. Another gener- 
alization is due to Eq. ((TO)) . according to which the same 
dependence for Lioc is expecting to occur when the dif- 
ference A = Uadb — \nbdb\ is sufficiently small, A <C cr^. 

Our results contribute to the theory of bilayered struc- 
tures with an inclusion of left-handed materials, and may 
be used in experimental realizations of structures with 
specific properties of transmission. 
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